#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use db_parameters;
use ortho_species;
use ensembl_parameters;
use chromosomes;
use min_max;
use die_error;
use run_cmd;
use DBI;
use align_gaps;
use Data::Dump qw(dump);

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_prot_position
	

	Description:
		Obtain the positions of each protein relative to other transcript of the
		same gene
	Populates DB tables:
		prot_positions

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

use constant GENEID =>0;
use constant CHRMSM =>1;
use constant START_ =>2;
use constant END___ =>3;
use constant STRAND =>4;
use constant CENTRE =>5;
use constant ORDER_ =>5;


#_________________________________________________________________________________________

#	save_prot_positions

#		
#		saves the gene position per protein

#_________________________________________________________________________________________
sub save_prot_positions($)
{
	my ($dbh) = @_;
	print STDERR "\t\tSaving protein positions per gene to panda...\n";

	#
	#	get number of transcripts per line
	#
	my $data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
    SELECT 
            prot_id,
            cnt_transcripts
        FROM
            (
            SELECT
                gene_id,
                count(*) AS cnt_transcripts
            FROM 
                ens_id
            GROUP BY 
                gene_id
            ) AS a NATURAL JOIN 
            ens_id;
PL/SQLCMD
	my %transcripts_per_gene;
	$transcripts_per_gene{$_->[0]} = $_->[1] for (@$data);


	#
	#	get transcripts in order of position
	#
	$data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
    SELECT 
            prot_id, 
            gene_position 
        FROM 
            gene_loci NATURAL JOIN 
            ens_id LEFT JOIN 
            prot_coding_structures pcs USING (prot_id)  
        ORDER BY 
            species, 
            gene_position, 
            (pcs.finish + pcs.start) / 2 ;
PL/SQLCMD
	$dbh->do("TRUNCATE taxon.prot_positions") or die;
	$dbh->do("VACUUM ANALYSE taxon.prot_positions") or die;
	$dbh->do("COPY taxon.prot_positions (prot_id, prot_position) FROM STDIN") or die;

	my $curr_gene_position = 10000000000000;
	my $sub_position = 0;
	for my $line(@$data)
	{
		my ($prot_id, $gene_position) = @$line;
		if ($gene_position != $curr_gene_position)
		{
			$sub_position = 0;
			$curr_gene_position = $gene_position;
		}
		my $digits = $transcripts_per_gene{$prot_id} > 1 ? 
				(int(log ($transcripts_per_gene{$prot_id} - 1) / log(10)) + 1) : 1;
		$gene_position += sprintf('0.%0'.$digits.'s', $sub_position++);
		$dbh->pg_putline($prot_id. "\t". $gene_position. "\n");
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;
	$dbh->do("VACUUM ANALYSE taxon.prot_positions") or die;
}















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# Main logic

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;

save_prot_positions($dbh);


print STDERR "\tCompleted\n";


$dbh->disconnect();


